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We provide molecular dynamics simulation of heat transport and thermal energy diffusion in one-dimensional 
molecular chains with different interparticle pair potentials at zero and non-zero temperature. We model the 
thermal conductivity (TC) and energy diffusion in the coupled rotator chain and in the Lennard-Jones chain 
either without or with the conhning parabolic interatomic potential. The considered chains without the conhning 
potential have normal TC and energy diffusion, while the corresponding chains with the confining potential are 
characterized by anomalous (diverging with the system length) TC and superdiffusion of energy. We confirm 
in such a way that, surprisingly, the confinement makes both heat transport and energy diffusion anomalous 
in low-dimensional phononic systems. We show that the chain, which has a finite TC, is also characterized 
by the normal energy diffusion in the thermalized chain (at non-zero temperature), while the superdiffusion of 
thermal energy occurs in the thermalized chains with only anomalous TC. We present the arguments, supported 
by our simulations, that the scaling relation between the exponents in time dependence of the mean square 
displacement of thermal energy distribution and in length dependence of the anomalous TC is not universal 
and can be different, depending on the main mechanism of energy transport: by weakly-scattered waves or by 
noninteracting colliding particles performing Levy flights. 


Phonon thermal conductivity in low-dimensional systems 
and at nanoscale presents a challenge both from theoretical 
and experimental points of view [1,2]. One of important prob¬ 
lems in this field is the origin of anomalous, diverging with the 
system size, thermal conductivity (TC) of one-dimensional 
(ID) lattice models [3-10] and quasi-lD nanostructures and 
polymers [11-13]. The length-dependent TC was experimen¬ 
tally observed in carbon nanotubes [14] and in suspended 
single-layer graphene [15]. The general belief is that TC is 
anomalous in ID momentum-conserving systems [16]. Nev¬ 
ertheless, a chain of coupled classical rotators presents an ex¬ 
ample of translationally-invariant (isolated) and correspond¬ 
ingly momentum-conserving ID periodic system with finite 
TC [17, 18]. In a recent paper [19], normal heat transport 
was reported for ID momentum-conserving systems with the 
Lennard-Jones, Morse, and Coulomb potential. It was shown 
in [19] that the convergence of TC is provided by phonon 
scattering on the locally strongly stretched loose interatomic 
bonds at low temperature and by the many-particle scattering 
at high temperature. 

While the research on TC in low-dimensional systems is 
very intensive, the peculiarities of energy diffusion have been 
much less studied. Usually the nonthermalized chain has been 
considered with only small central portion being thermalized, 
and thermal energy propagation in the chain was studied af¬ 
terwards. Such modeling can lead to anomalous (super-) dif¬ 
fusion of the thermal energy, and sometimes it is concluded 
that such ID system possesses anomalous thermal conductiv¬ 
ity [20] . For example, superdiffusion of energy was obtained 
in Ref. [20] in the chain with the Lennard-Jones (LJ) inter¬ 
particle potential, in which the normal heat transport has been 
revealed previously [19, 21]. This suggests the incorrectness 
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of the simulations of energy diffusion because the latter was 
simulated in Ref. [20] in the chain at zero background temper¬ 
ature. Any anharmonic chain at low temperature behaves as 
an almost harmonic chain, in which the heat transport is car¬ 
ried out by ballistic phonons and TC diverges with the chain 
length. Therefore it is not surprising that in the nonthermal¬ 
ized chain (at zero background temperature) the energy prop¬ 
agates ballistically. But this does not imply that the superdif¬ 
fusion of energy will also hold in the thermalized chain (at 
non-zero background temperature). 

On the other hand, the chains with a (parabolic or quar- 
tic) confining pair potential, which does not allow for bond 
dissociation, are expected to possess anomalous thermal con¬ 
ductivity, diverging with the chain length [19]. The celebrated 
Fermi-Pasta-Ulam potential belongs to such type of the pair 
potentials. In this paper, we model both the TC and diffusion 
of energy in anharmonic chains at zero and non-zero temper¬ 
ature. We show that the chain, which has a finite TC, is also 
characterized by the normal diffusion of energy in the ther¬ 
malized chain, while the anomalous superdiffusion of ther¬ 
mal energy occurs in the thermalized chains with only anoma¬ 
lous (diverging with the system length) TC. To show this, we 
model TC and energy diffusion in the coupled rotator chain 
(a chain with periodic interatomic potential) and in LJ chain 
either without or with the parabolic confining pair potential. 
Such parabolic confining pair potential is introduced, for in¬ 
stance, on top of the shallow LJ potential to hold together the 
constituent cellular particles within a cell [22]. The consid¬ 
ered chains without the confining potential are characterized 
by the normal diffusion of energy in the thermalized chain, 
while the corresponding chains with the confining pair poten¬ 
tial are characterized by the superdiffusion of energy. We con¬ 
firm in such a way that, surprisingly, the confinement makes 
both heat transport and energy diffusion anomalous in low¬ 
dimensional systems. We present the arguments, supported 
by our simulations, that the scaling relation between the ex¬ 
ponents in time dependence of the mean square displacement 
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of thermal energy distribution and in length dependence of the 
TC, which is discussed in connection with anomalous heat 
transport and superdiffusion of energy in low-dimensional 
systems, see, e.g.. Refs. [23-27], is not universal and can be 
different, depending on the main mechanism of energy trans¬ 
port - either by weakly-scattered waves or by noninteracting 
colliding particles performing Levy flights. 

We consider the periodic molecular chain consisting of 
N = L/a unit cells, where L is the chain length, a is a lattice 
period. In a dimensionless form, the Hamiltonian of the chain 
can be written as 

N ^ N-1 

H = ^ ^ 2^” ^ ^ V{Xn+l ~ Xn), (1) 

n—0 n—0 

where Xn is the displacement on the nth particle from its equi¬ 
librium position at na, V(r) is a dimensionless pair interac¬ 
tion potential between nearest neighbors normalized by the 
conditions L(0) =0 and L'(0) = 0- 

To simulate the heat transfer in the chain, we use the 
stochastic Langevin thermostat. We consider a finite chain 
of + N + iV_ unit cells, and take the chain with fixed 
ends. We put the N- left boundary particles in the Langevin 
thermostat with temperature r_, and particles in the right- 
hand edge Langevin thermostat with a temperature of r+. The 
corresponding system of equations of motion of the chain is: 

Xn = -dff/dxn - JXn + n<N+, 

Xn = -dH/dxn, n = N+ + l,...N+ + N, (2) 

Xn = -dH/dXn-lXn + in^ n> N++ N, 

where 7 is relaxation coefficient of the particle velocity, 
are random forces which simulate the interaction with the 
white-noise thermostat normalized by the conditions 

W) = 0, = 2jT±6nkSih - h). 

We integrate numerically the Langevin equations of mo¬ 
tion (2) by employing the Verlet velocity method with the 
step At = 0.02. After some integration time to (this value 
depends on the chain length between the thermostats), we ob¬ 
serve the formation of a temperature gradient and constat heat 
energy flux in the central part of the chain. After the station¬ 
ary heat flux is established, we can find the temperature dis¬ 
tribution by using the relations for r„ = {Xn)t and station¬ 
ary heat flow along the chain = —{xnV'{xn — Xn-i))t- 
The following values were used in the numerical simulation: 

r± = (1 ± 0.1)r, 7 = 0.1, N± = 40, A = 20, 40, 80, 
..., 20480. The detailed description of the Langevin equations 
and justification of the method are presented in [19]. 

In the steady-state regime, the heat flux through each cell in 
the central part of the chain should be the same, i.e., J„ = J, 
A_ + 1 < n < A_ -f A -f 1. Almost linear gradient of 
temperature distribution is established in the central part of 
the chain, so we can define TC as 

k{N) = J{N - 1)/{Tn^+i - Tn^+n). (3) 

In our modeling of the TC, we use the following pair inter¬ 
action potentials: 

L(r) = 1 — cos(r), a = 27r, (4) 



Figure 1: Thermal conductivity k versus dimensionless chain length 
N = L/a for the chain with the periodic potential (4) at normalized 
temperature T = 0.15, 0.2, 0.3 (curves 1, 2, 3), for the chain with 
the combined potential (5) at T = 0.3 (curve 4), for the chain at 
T = 0.002 with the Lennard-Jones potential (6) (curve 5, which 
shows for the convenience 0.15k) or with the combined potential (7) 
(curve 6). 

- the periodic potential, 

L(r) = 1 — cos(r)-I-0.25r^, a = 27r, (5) 

- the sum of the periodic and parabolic confining potentials, 

V{r) = Ae\{a/{l+r)f-l/2 ]\a = 1, a = e = 1/72, 

( 6 ) 

- the LJ potential, 

V{r) = 4e[(cr/(l -f r)f - l/2]2 -p 0.25r'^, (7) 

- the sum of the LJ and parabolic confining potentials. Note 
that one has V"{Q) = 1 for potentials (4) and (6) and therefore 
speed of sound (small-amplitude phonons) in the chain is Vg = 
a, while for potentials (5) and (7) one has L"(0) = 1.5 and 
therefore Vg = aVl-5 = 1.2247a. 

Dependence of TC on dimensionless chain length k(A) is 
presented in Fig. 1. As one can see from this figure, k(A) in 
the chain with the periodic potential (4) depends on tempera¬ 
ture: the convergence of k(A) for A 00 is slower for lower 
temperature. Thermal conductivity saturates for A = 2560 at 
T = 0.3, for A = 20480 at T = 0.2, and there is no saturation 
of the k(A) at the maximal used A = 20480 at T = 0.1 (and 
one needs A ~ 10^ to reach the saturation at such tempera¬ 
ture). This feature is related with the increase of phonon mean 
free path with the decrease of temperature because the satura¬ 
tion of k(A) occurs only when the chain length L = Na 
reaches or exceeds phonon mean free path. 

The chain with the combined interatomic potentials (5) and 
(7) is characterized by anomalous heat transport. Here TC 
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monotonously diverges with N as k{N) ^ N°‘. For the po¬ 
tential (5) at r = 0.3, the exponent a = 0.705, for the poten¬ 
tial (7) at T = 0.002, the exponent a = 0.76. TC of the chain 
with the LJ potential ( 6 ) saturates for N ~ 10^ at T = 0.002, 
as one can see in Fig. 1. 

Now we analyze how the k( 7V) dependence is related with 
thermal energy propagation at finite temperature (Tq > 0). To 
this end, we consider the chain with N atoms with the fixed 
ends. We thermalize the chain with the use of the Langevin 
heat bath such that the first Ni atoms have the temperature 
Ti, while the rest of the atoms have the temperature Tq < Ti. 
We consider the Langevin equations of motion with random 
forces; 

Xn=-dH/dXn-JXn+^n, 0 < U < N, ( 8 ) 

when lo = 0 , xn = 0 , and the random forces are 6- 
correlated as 

(Cn(fi)^fc(f2)) = 27T„(5„fc(5(f2 - h), 

with Tn = Ti for n = 1,..., Ni, and = Tq for n = Ni + 
I, ...,N — 1. Now we take the (classical) zero-energy initial 
condition {a;„( 0 ) = 0 , XniO) = and integrate the 

Langevin equations ( 8 ) from f = 0 to L = 2 O /7 = 200. In 
result we get a random realization of the initial state of the 
thermalized chain {xni.ti),Xn{ti)}^^-l with T = Ti in its 
left end and t = Tq < Ti in the rest of the chain. 

Then we model the propagation of the thermal energy in the 
chain, which starts from its left end with higher temperature. 
To perform this, we integrate equations of motion of the atoms 
without their interaction with the thermostat, 

Xn = —dH/dxn, 0 < n < N, (9) 

when a;o = 0 , xat = 0 , with the initial condition which were 
obtained after the integration of the Langevin equations of 
motion ( 8 ). In order to analyze the propagation of the ther¬ 
mal energy during the long time in the chain with N = 20000 
atoms, we take the number of the end atoms = 40 <C N. 
We trace the time dependence of the distribution along the 
chain of the temperature T„(f) = (i^(f)) and of the energy 
En{t) = {Xn{t)/2 + V{xnit) - Xn-i{t))), where the aver¬ 
aging is taken over the independent realizations of the initial 
thermalized state of the chain. Time evolution of the tem¬ 
perature distribution in the chain with the periodic potential 
is shown in Fig. 2. It is worth noting that for the analysis 
of thermal energy distribution in the thermalized chain, with 
Tq > 0, one needs to perform the averaging over a large num¬ 
ber 10 ®) of independent realizations of the initial thermal¬ 
ized state [while for the zero-energy initial state, with Tq = 0, 
the averaging can be performed over significantly lower num¬ 
ber 10 ®) of independent realizations of the initial state]. 

We consider first the chain with the periodic interatomic po¬ 
tential (4). The left end of the chain is thermalized at Ti = 1, 
and the rest of the chain - at Tq, when 0 < Tg <Ti. As one 
can see in Fig. 2, the thermal energy distribution substantially 
depends on the chain temperature Tq. For Tq = 0, the thermal 
energy propagates ballistically, with the sound speed, along 



Figure 2: Propagation of temperature distribution Tn in a chain with 
the periodic potential (4) for Ti = 1 and Tq = 0 (a). To = 0.15 (b). 
To = 0.2 (c). Temperature distributions are shown for characteristic 
delay time dt = 1662.5, dashed lines show sound cone n = Vst/a. 


the chain. While for Tq = 0.15, only relatively small part of 
thermal energy propagates ballistically, which then dissipates 
in the chain. The main part of thermal energy is concentrated 
at the left, high-temperature, end of the chain and spreads 
slowly towards the right end. For the higher temperature of 
the chain, Tq = 0 . 2 , there is no ballistic energy propagation 
along the thermalized chain: the energy distribution spreads 
diffusively from the very beginning. 

The unilateral spreading can be described by the follow¬ 
ing mean square displacement of thermal energy distribution 
(MSDTED) of the excess energy, initially placed in the site 
n = 1 of the chain with uq = 0 , un = 0 : 

Af-l 

(Ax^)(f) = ^ (n - l)^e„(f). (10) 

n—1 

Heree„(f) = (£'„(f)—i?o)/£’is a discrete distribution of nor¬ 
malized excess energy in the chain, E = — Eq) 

is a constant total excess energy, Eq = Tq is the average par¬ 
ticle energy in the thermalized chain at temperature Tq, see 
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Figure 3: Time dependence of MSDTED, Eq. (10), in the chain 
with the periodic potential (4) for the chain temperature To — 0, 
0.15, 0.2, and 0.3 (curves 1, 2, 3, and 4). Dotted lines show the linear 
dependencies. 


also Refs. [27, 28]. Mean square displacement can also be 
introduced for the discrete temperature distribution along the 
chain, when T„(f) = and both mean square displace¬ 

ments, of thermal energy and temperature distributions, have 
the same time dependence. 

Time dependence of MSDTED in the chain with the peri¬ 
odic potential (4) is shown in Fig. 3. As one can see in this fig¬ 
ure, at zero and low chain temperature, Tq = 0 and Tq = 0.15, 
MSDTED grows quadratically with time, (Ax^) oc which 
corresponds to the ballistic energy propagation. At higher 
temperature Tq = 0.2, during the short delay time t < 10^ 
MSDTED also grows quadratically with time, but for longer 
delay t > 10"*^ MSDTED changes to become linear in time, 
which corresponds to the normal energy diffusion. As one can 
also see in Fig. 3, the normal energy diffusion starts earlier for 
the higher temperature of the chain which is a consequence of 
the enhanced phonon scattering by thermally-activated anhar- 
monicity of the chain. From Fig. 3 follows that the used chain 
length N is not enough long to model the transition from the 
ballistic to normal diffusion regime of energy propagation at 
temperature Tq = 0.15. We also note that this length of the 
chain is not enough to model the normal TC in it at this tem¬ 
perature as well. We can conclude from the comparison of 
Figs. 1 and 3 that the minimal chain length Nmin, at which the 
normal TC is established, is directly related with the minimal 
delay time which is needed for the normal energy diffu¬ 
sion to be established in the same chain; ~ 

Similar time dependence of MSDTED is revealed in the 
chain with the LJ interatomic potential, see Figs. 4 and 5. We 



Eigure 4: Propagation of energy En in the chain with the Lennard- 
Jones potential (6) for Ti = 0.01, To = 0.0 (a) and To — 0.002 (h). 
Energy distributions are shown for characteristic delay time dt = 
1000, dashed lines show sound cone n = Vatja. 


consider the chain with the same N = 20000, when the first 
A) = 40 atoms in its left end are thermalized at T/ = 0.01. 
As one can see in Fig. 4 (a), for Tq = 0 the energy spreading 
can be separated into the two parts, which propagate with the 
sound and supersonic speeds, respectively. The first part cor¬ 
responds to energy transfer by small-amplitude linear lattice 
waves (phonons) while the second one corresponds to energy 
transfer by supersonic kinks (compression acoustic solitons) 
produced by hard compression component of the LJ inter¬ 
atomic potential [29-31]. For Tq = 0 MSDTED increases 
quadratically with time, (Aa:^) ^ see line 1 in Fig. 5, 
which corresponds to the ballistic energy propagation in the 
zero-temperature (nonthermalized) chain. 

The picture of energy propagation changes in the thermal¬ 
ized chain. For Tq = 0.002, the part of the initial energy, 
which propagates with supersonic speed, starts to dissipate in 
the chain (because the acoustic solitons have finite mean free 
path in the thermalized LJ chain), and the normal diffusion of 
thermal energy is established in the chain, see Fig. 4 (b). Time 
dependence of MSDTED shows that for t > 5000 the ini¬ 
tial ballistic energy propagation (superdiffusion) is replaced 
by the normal energy diffusion when (Ax^) grows linearly 
with time, see line 2 in Fig. 5. In the chain with the combined 
potential (5) or (7), the anomalous superdiffusion of energy 
takes place and MSDTED grows as a power function of time: 
(Ax^) oc with the exponent 1 < /3 < 2. In the chain with 
potential (7), /3 = 1.61 at Tq = 0.002, see line 3 in Fig. 5; in 
the chain with potential (5), /3 = 1.544 at Tq = 0.3, see line 2 
in Fig. 6 . Again the minimal delay time, at which the normal 
energy diffusion is established, is consistent with the minimal 
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Figure 5; Time dependence of MSDTED, Eq. (10), in the chain with 
the Lennard-Jones potential (6) with temperature Ti = 0.01 of the 
left end and temperature To = 0 and To = 0.002 (solid lines 1 and 2) 
of the rest of the chain, and in the chain with the combined potential 
(7) for To — 0.002, Ti = 0.01 (solid line 3). Dashed lines 1, 2 and 
3 give the quadratic 1.41f^, linear 6800(f — 2000) and power-law 
31(f — 300)^ ®^ dependencies. 


chain length, at which the normal TC is established, cf. Figs. 
5 and 6 with Fig. 1 . 

The value of the exponent /3 can be related with the ex¬ 
ponent a in the length dependence of anomalous TC, when 
K oc N°‘. By the definition of the thermal energy diffusion 
De and TC k = cDe coefficients, where c is the specific heat 
density, MSDTED grows as a power function of time as 

(Aa;^) oc DE{N)t oc N°‘t, (11) 

(Ax^) oc (12) 

where De{N) = k{N)/c oc iV“. Below we show that (at 
least) two different relations between a and (3 can be obtained 
from Eqs. (11) and (12) under different assumptions. Accord¬ 
ing to its definition (10), (Ax^) scales as the square of the 
dimensionless length. Then under the assumption 

(Ax^) oc 7V^ (13) 

we get from Eqs. (11) and (12) that 

a = 2-2//3. (14) 

Relation (13) corresponds to the assumption that the regime 
of the thermal diffusion is reached only in the chain, whose 
length reaches (or exceeds) the effective phonon mean free 
path. But under the different assumption 


(15) 



Figure 6; Time dependence of MSDTED, Eq. (10), in the chain with 
the combined potential (5) with temperature T; = 1 of the left end 
and temperature To — 0 and To = 0.3 (solid lines 1 and 2) of the rest 
of the chain. Dotted lines 1 and 2 give the dependencies 0.7t^ and 
12.5(f — 110)^ ®"'"'. Solid lines 3 and 4 give the dependencies for the 
chain with the periodic interatomic potential (4) for Tq = 0.2 and 
To = 0.3. Dashed lines 3 and 4 give the dependencies 230(f — 800) 
and 45 (f - 964). 


we get from Eqs. (11) and (12) that 

a = /3-l. (16) 

Relation (15) corresponds to the assumption that the energy 
carriers launched from the hotter side of the chain will (ballis- 
tically) reach the colder side and reflect back beyond the delay 
time t. Both scaling relations (14) and (16) imply that normal 
energy diffusion {f) = 1) leads to the normal (non-divergent) 
TC (a = 0), while the superdiffusion (/3 > 1) corresponds to 
anomalous (a > 0) TC. Relation (14) was suggested in Refs. 
[8, 23], while the relation (16) was obtained in Refs. [24-27] 
in the specific case of billiard-like ID models in which non¬ 
interacting particles undergo Levy flights. It is worth noting 
that the assumption (13) was implicitly used in derivation of 
Eq. (14) in Ref. [8], see also Ref. [32], and the assumption 
(15) was explicitly used in derivation of Eq. (16) in Ref. [27]. 
Moreover, the scaling relation (14) can be derived from the 
analysis of the diffusion with a position-dependent diffusion 
coefficient, see, e.g.. Ref. [33], under the same assumption 
(13). 

The main conclusion of our analysis is that neither of the 
relations (14) and (16) is the universal relation, which can 
be applied to all the nonlinear systems with anomalous heat 
transport and superdiffusion of thermal energy. We can com¬ 
pare this conclusion with the known conjecture on the absence 
of the unique velocity-correlation function in turbulent flow, 
which is universal for all relevant scales and types of flow, 
see Refs. [34, 35]. Our modeling of anomalous thermal con¬ 
ductivity and superdiffusion of thermal energy in the chains 
with the combined interatomic potentials (5) and (7) confirms 
with high accuracy the relation (14): we have a = 0.705, 
P = 1.544 for the chain with the potential (5), and a = 0.76, 


t oc Na/vs, 
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P = 1.61 for the chain with the potential (7). On the other 
hand, the relation (16) was confirmed in the study of heat 
transport in the billiard-like ID system containing colliding 
particles with two different masses at Tq = 0 [25, 36]. We 
relate the difference between Eqs. (14) and (16) with the fact 
that Eq. (14) is applied mostly to the ID lattices of the cou¬ 
pled (anharmonic) oscillators, in which heat is transported by 
weakly-scattered waves (phonons), while Eq. (16) is applied 
mostly to the billiard-like ID systems, in which heat is trans¬ 
ported by noninteracting parficles performing Levy flights. 

In conclusion, we show that the normal thermal conduc¬ 
tivity is always accompanied by the normal energy diffusion 
in the thermalized anharmonic chains, while the superdiffu¬ 
sion of energy is inherent in the thermalized chains with only 


anomalous heat transport. We confirm that the confining in¬ 
terparticle potential makes both heat transport and energy dif¬ 
fusion anomalous in low-dimensional phononic systems. We 
show that the scaling relation between the exponents in time 
dependence of the mean square displacement of thermal en¬ 
ergy distribution and in length dependence of anomalous ther¬ 
mal conductivity is not universal and can be different, depend¬ 
ing on the main mechanism of energy transport: either by 
weakly-scattered waves or by noninteracting colliding parti¬ 
cles performing Levy flights. 

The authors are grateful to the Joint Supercomputer Center 
of the Russian Academy of Sciences for the use of computer 
facilities. 
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